Simulation of a Single Polymer Chain in Solution by Combining Lattice Boltzmann 

and Molecular Dynamics 



«4H 

o 



Patrick Ahlrichs, Burkhard Diinweg 
Max Planck Institute for Polymer Research, Ackermannweg 10, D-55128 Mainz, Germany 

(February 1, 2008) 

In this paper we establish a new efficient method for simulating polymer-solvent systems which 
combines a lattice Boltzmann approach for the fluid with a continuum molecular dynamics (MD) 
model for the polymer chain. The two parts are coupled by a simple dissipative force while the 
system is driven by stochastic forces added to both the fluid and the polymer. Extensive tests of 
the new method for the case of a single polymer chain in a solvent are performed. The dynamic and 
static scaling properties predicted by analytical theory are validated. In this context, the influence 
of the finite size of the simulation box is discussed. While usually the finite size corrections scale 
as (L denoting the linear dimension of the box), the decay rate of the Rouse modes is only 
subject to an finite size effect. Furthermore, the mapping to an existing MD simulation of the 
same system is done so that all physical input values for the new method can be derived from pure 
MD simulation. Both methods can thus be compared quantitatively, showing that the new method 
allows for much larger time steps. Comparison of the results for both methods indicates systematic 
deviations due to non-perfect match of the static chain conformations. 

PACS numbers; 02.70.Ns, 05.10.Gg, 47.11. +j, 83.10.Nn 
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I. INTRODUCTION 



The complexity and variety of soft condensed matter 
is largely due to the fact that IpjPgth scales of different 
orders of magnitude are presentuD. When dealing with 
polymers in computer simulations, one therefore often in- 
tends to analyze the scaling behavior, where the natu 
of the underlying chemistry becomes unimportanlld' 
When constructing models for these systems it is crucial 
to coarse-grain the details and to keep the relevant length 
scales in order to observe the phenomena one is interested 
in. Since bead-spring models in MD simulations are an 
appropriate means to yield the right universal laws, they 
have been widely used to simulate the scaling behavior of 
polymera_and much progress has been made using these 
modelsOLl. 

While in some systems, e. g. in highly concentrated 
solutions or in melts, the dynamic properties are not af- 
fected by the solvent — such that these can be simulated 
by conventional bead-spring models without explicitly 
taking into account the solvent — there are many phe- 
nomena in polymer science where the influence of the 
solvent on the polymer dynamics cannot be neglected. 
For example, in dilute or semi-dilute polymer solutions, 
the dynamical behavior is changed and even dominated 
by hydrodynamic interaction between different parts of 
the polymers. This eventually leads to a long-range in- 
teraction which is mediated by the solven#0. With this 
paper, we want to provide a new efficient method for 
the simulation of polymer systems where hydrodynamics 
plays a role. The idea is to focus on the really necessary 
parts only, i. e. the hydrodynamics of the solvent and the 
(Brownian) motion of the polymer chains, thereby try- 
ing to keep the computational costs at a minimum. Our 



test case is the dynamics of a single chain in a solvent. 
This problem haSp-continuously attracted the attention 
of MD resejaxchersQ^EI, mainly because existing analytical 
theoriesllEI E3 rely on uncontrollable assumptions that can 
be tested using computer simulations. 

Simulating such systems by MD is only possible if one 
introduces explicit solvent particles. Hence one has to 
face the problem that almost all CPU time goes into the 
propagation of the solvent particles, while one is mainly 
interested in the chain properties. However, there are 
also other computational methods than MD available for 
soft condensed matter systems where hydrodynamics is 
important, not only in the field of polymers but for ex- 
ample also in colloidal susppisians. These include Brow- 
nian Dynamics simulai;ipagl3il3, and Dissipative Parti- 
cle Dynamics (DPD)EZl"E^. Both of them have inherent 
strengths, but also some disadvantages: The first tech- 
nique must face the problem that the algorithm scales 
as the cube of the number of particles, and the latter 
(like MD) simulates the solvent particles explicitly, lead- 
ing to simulations of several thousand particles even for 
a single chain of, say, 30 monomers. Compared to MD, 
DPD has the advantage of much larger timeiieps, mainly 
because of the use of very soft potentialsllj. A lot of 
progress in th|e-|tlieoretical framework of the method has 
been achievedE£rE3, but some practical problems remain, 
like the time step-dependent temperature and the small 
Schmidt numberllSI. Receiiijly, however, some effort has 
been made to fill this gapO. |_. 

In this paper we use a recently proposed methodEj 
that couples a lattice Boltzmann approach for the fluid 
to bead-springj^olymer chains. The lattice Boltzmann 
method (LBM)c3'E3 was developed to simulate hydrody- 
namics on a grid. The LBM was shown to be an effective 
and fast method for simulating fluid flows, comparable to 
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finite-difFerencJni or spectral methodfl. rtjadd applied 
the LBM successfully to colloidal systemsE3o: The col- 
loidal particles are simulated as hard spheres by using 
stick boundary conditions. This leads to a very efficient 
algorithm: Its CPU cost scales linearly with the num- 
ber of particles, and it uses a "minimal" rnodel to sim- 
ulate the fluid. Besides, Ladd also showeda that fluc- 
tuations can be incorporated into the LBM inptiie spirit 
of Landau-Lifshitz fluctuating hydrodynamicscj, which 
is essential if one wants to investigate Brownian motion. 

Now one might think of a direct application of Ladd's 
method to polymer-solvent systems. However, using 
hard spheres to model the monomers is not necessary 
here, as rotational degrees of freedom as well as stick 
boundary conditions are not relevant: On the large 
length and time scales we are interested in, like the ra- 
dius of gyration and the Zimm time of the polymer, it 
is sufficient that hydrodynamic interaction has evolved. 
The "microscopic details" of the coupling should then 
not play a role. In this spirit, we couple the LBM to 
bead-spring polymer chains by a simple friction ansatz, 
thereby treating the monomers as point particles for the 
fluid. We will show that this ansatz is sufficient to sim- 
ulate both the static and dynamic scaling behavior of 
the polymer. The simulation of the fluid by LBM rather 
than explicit particles and the simple friction ansatz lead 
to a large speedup in computer time of about a factor of 
20 when compared to pure MD, or even more if one is 
willing to be satisfled with less accurate data. 

Additionally, we map our method to a pure MD simu- 
lation, i. e. we show how to determine all physical input 
values from the results of MD, allowing us to compare our 
results t»|an existing MD simulation with explicit solvent 
particlesEI. In other words, the fluid in the new method 
can be viewed as a coarse-grained MD fluid, and there 
exists a well-defined procedure for how to do the coarse- 
graining. Of course, in using such a mesoscopic approach 
it is no longer possible to include detailed chemistry like 
in atomistic MD simulations. This is, however, a quite 
common feature of mesoscopic simulation methods; DPD 
simulations do not include atomistic details either. 

The remainder of this article is organized as follows: 
We outline the method in Section ||, and present the nu- 
merical results in Section p^ , which arc compared to pure 
MD in Section 0. In Section wc conclude with some 
final remarks and an outlook to further studies. 



II. THE SIMULATION METHOD 



A. The Lattice Boltzmann Method for the Solvent 



The lattice Boltzmann method is a discrete formula- 
tion of the Boltzmann equation on a lattice, leading to 
the Navier-Stokes equations in the incompKessible limit 
by means of a Chapman-Enskog expansionOEa. It has 



been successfully applied to a variety of fluid flow prob- 
lems, and it is especially well-suited for complex fluids 
because of the possibility of straightforward implementa- 
tion of complex boundaries. The central quantity of the 
algorithm is ni{r,t), the number of particles in a volume 
at the grid point r at time t, which have the veloc- 
ity Ci^ (i = 1, ..,6), where a is the lattice spacing, r the 
time step and a vector leading to the ith neighbor on 
a grid with unit lattice constant. The evolution equation 
for ni(r,t) is the lattice Boltzmann equation 



ni{r + c^a, t + t) = nj(r, t) 



(1) 



The last term expresses the relaxation of towards a 
local pseudo-equilibrium, which resembdes a Bhatnagar- 
Gross-Krook (BGK) collision operatoicll in the contin- 
uum Boltzmann equation. The constant matrix Ly can 
be interpreted as the scattering between particle popu- 
lation i and j. Its eigenvalues can be determined from 
physical and numerical arguments, such that its |explicit 
form is not necessary for the simulation algorithmc£l. The 
local pseudo-equilibrium distribution n^'(p, u) depends 
on the density p(r, t) = n.i(r, t)^/a'^ and fluid current 
j(r, i) = pu = X). ni(r,t)cjp/(ra2) only. Here, fi is the 
mass of a fluid particJe. The usual functional form for 
71^'^ (p, u) is assumedO: 



(p, U) = p{Ay+ Bq {Ci ■ U) + CgU^ + Dq (c, • u)^ 



(2) 

The coefficients Aq, Bq, Cq and Dq (which depend on 
the sublattice q, i. e. the magnitude of c^) are deter- 
mined to reproduce the correct macroscopic hydrody- 
namic behavior. Note that this is contrary to contin- 
uum kinetic theory, where the Maxwell-Boltzmann dis- 
tribution is determined by entropy considerations and 
the Navier-Stokes equations_.fQLllow naturally by the 
Chapman-Enskog expansiorSH. Hence it is called 
pseudo-equilibrium. Explicit values for the coefficients 
Aq,Bq,Cq and Dq are known for different latticeg^. 

Here, we implement the 18-velocity model of Ref. |29| , 
which corresponds to the D3Q18 model in the nomencla- 
ture of Ref. |3^. The set of consists of the 6 nearest and 
12 next-nearest neighbors on a simple cubic lattice. Via 
a Chapman-Enskog expansion one can show that this 
model leads to the Navier-Stokes equations in the limit 
of small Knudsen and Mach numbersE3, and derive a re- 
lation between the kinematic viscosity v and the non- 
trivial eigenvalue A of Lij belonging to the eigenvector 
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T 



(3) 



In this paper, we always deal with low Reynolds num- 
ber flow, hence the linearized Navier-Stokes equations 
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are sufficient. For this reason, we neglect the nonlinear 
term in the equilibrium distribution i. e. we effec- 
tively set Cq =rQq = 0, thus obtaining a simpler and 
faster algorithmcj. 

Fluctuations-can be incorporated into the lattice Boltz- 
mann methocO. The central idea is to add fluctuations 
to the fluxes of the conserved variables, i. e. the stress 
tensor, and not to the hydrodynamic fields p and j. In 
this way, local mass and momentum conservation can be 
guaranteedEj. The fluctuating lattice Boltzmann equa- 
tion reads 



the LBM in dimensional form in the last section, rather 
than using the usual dimensionless lattice units. The 
equations of motion resulting from these potentials are 
integrated using the velocity Verlet algorithmE3 with a 
time step At. Note that there is a priori no need to set 
At = T and we will exploit this fact below. 

The polymer model has been lanplied successfully to 
the simulation of many, systemsoU including a single 
chain in explicit solventQ, so that we can compare chain 
properties in using these potentials. 



ni(r -I- Cifl, t + t) = nj(r, t) 

b 

+ {nj{r,t) -n''/{p,u)) +n'^{r,t) 



(4) 



with the stochastic term 
<ir,t) 



a/3 



(5) 



The random stress fluctuations a'^^ are assumed to have 



white noise behavior 
{a',f,{r,ty^s{r',t'))^A6,,,6u' 



(6) 



By solving the resulting discrete Langevin equation 
for the_current one finds the fluctuation-dissipation 
relatiorO for this system; the noise strength A is given 
by 



A 



2r]kBTX^ 



(7) 



where rj = vp is the dynamic viscosity. 

The LBM was tested extensively, compared to other 
Navier-Stokes solvers and found to have comparable 
speed and accuracy (see for example Refs. |2^,| 



B. The Bead— Spring Model for the Polymer Chain 

The polymer model consists of repulsive Lennard- 
Jones monomers |j:onnected via non-harmonic springs 
(FENE potential)^: 



14j = 4e 
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FENE 
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(r<2i/V) (8) 
(r < i?o). 



In order to model the excluded volume effect the 
Lennard-Jones potential acts between all monomers. As 
usual, the parameters e, a and the mass m of the 
monomer define our unit system. Therefore we wrote 



C. Coupling of Fluid and Monomer 



As mentioned above, for the length and time scales of 
the polymer chain, the "microscopic" details of the cou- 
pling should not play a role, as long as one assures that 
hydrodynamics evolves in the fluid on time scales faster 
than the diffusion time scale of the monomers. It is not 
necessary to resolve the shape of the monomer for the 
fluid. Thus, we can treat one monomer as a point par- 
ticle. In analogy to the Stokes formula for a sphere in a 
viscous fluid, we assume the force on the monomer ex- 
erted by the fluid to be proportional to the difference of 
the velocity of the monomer V and the fluid velocity u 
at the monomer's position. 



-C[V-u(R,i)], 



(9) 



Here, C is a proportionality coefficient which we will refer 
to as the "bare" friction coefficient. This ansajtp has also 
been used in the simulation of sedimentationtj. 



a- Ay 



Ay 



(0,1) 


(1,1) 


(0,0) 


(1,0) 



Ax 



a- Ax 



FIG. 1. Illustration of the quantities used for the coupling 
of monomer and fluid (in two dimensions) . The figure shows a 
sketch of a monomer surrounded by the elementary cell of the 
four nearest neighbor grid points, a is the lattice constant. 
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Because the fluid velocity is only calculated at the dis- 
crete lattice sites in the simulation, one has to interpolate 
to get u(R, t) at the monomer's position. We implement 
a simple linear interpolation using the grid points on the 
elementary lattice cell containing the monomer: Denot- 
ing the relative position of the monomer in this cell by 
(Ax, Ay, Az), with the origin being at the lower left front 
edge (see Fig. |l|), we can define 

^(0,0,0) = (1 - Aa;/a)(l - Ay/a)(l - Az/a), (10) 
^(1,0,0) = Ax/a • (1 - Ay/a){l - Az/a), 

etc. The formula for the linear interpolation then reads 
u(R,<)= ^5ru(r,t) (11) 

rSng 

where ng denotes the grid points on the considered ele- 
mentary lattice cell. 

In order to conserve the total momentum of fluid and 
monomer we have to assign the opposite force to the fluid 
in that cell. Note that then the interaction is purely local. 
In particular, the force density —Ffi/a^ which is to be 
given to the fluid leads to a momentum density transfer 
per MD time step At of 

'i,rGng 

The last equation has to be satisfied for the change in 
the number of particles Arii of the grid points on the ele- 
mentary lattice cell in order to exchange the momentum 
density Aj. Besides, one must also ensure mass conser- 
vation in the fluid, 

J2 An,(r,t) = 0. (13) 

i,r6ng 

The way how to calculate the corresponding Arii at the 
nearest grid points is not unique; one possibility was pre- 
sented in Ref. HJ. Here, we follow a different approach 
which seems slightly more natural: For given hydrody- 
namic fields p{r,t) and j(r,t) at a certain grid point r, 
the equilibrium distribution can be calculated according 
to Eq. ||. The change in the equilibrium distribution at 
the points r £ ng due to the presence of the monomer can 
therefore be determined: p stays constant (mass conser- 
vation), while j ^ j -I- ^rAj. Here 6r is the fraction ( |l0| ) 
of the total Aj which is given to the specific grid point r. 
Therefore, by requiring that rii — n^"^ remains unchanged, 
we obtain 

An,{r) ^ BgSrAi ■ c„ (14) 

where again the nonlinear part of Eq. ^ has been ne- 
glected, consistent with our overall procedure. More 
accurate algorithms (which would however be computa- 
tionally more expensive) could be constructed, using the 



method proposed in Ref. however, this is not neces- 
sary for our purposes: Our simple approach is consistent 
with locality of the interaction, plus momentum conser- 
vation, and should therefore suffice to build up hydrody- 
namic interactions in the correct manner. 

As we discussed in Ref. |2j, one has to take care when 
adding stochastic terms to the system. Due to the dis- 
sipative nature of the coupling, it is necessary to incor- 
porate fluctuations to both the fluid and the monomers, 
i. e. to the LBM like in Eq. |^, and to the monomers by 
extending Eq. ^ to 

F/;--C[V-u(R,i)]+f. (15) 

Here f is a stochastic force of zero mean and 

{Ut)Uit')) = S{t - t')25^pkBT(:. (16) 

The momentum transfer to the fluid for the fluctuating 
case is calculated in the same way as described above 
without the fluctuations. For this reason, the total mo- 
mentum of fluid and polymer is conserved locally also in 
the fluctuating case. One can show analytically that with 
this method the fluctuation-dissipation relation holds for 
the continuum limit of the model, where the coupling to 
the LBM fluid is replaced by the analogous coupling to a 
Navier-Stokes fluid with thermal fluctuations of the flow 
field. For the velocities of the monomers and the fluid 
flow velocity, the equilibrium distribution is then given by 
the Maxwell-Boltzmann distribution, while the confor- 
mational statistics of the chain is given by the Boltzmann 
distribution, i. e. governed by the intra-chain potentials 
Vlj and Vfene, see Eq. |[ This should be contrasted 
with the MD case, where the potential due to the solvent 
particles has an additional influence. For the discrete 
case, one can check the fluctuation-dissipation relation 
by investigating the velocity relaxation of one (initially 
kicked) monomer in the fluid on the one hand, and the 
velocity autocorrelation, if fluctuations are added, on the 
other hand. The two quantities coincide for our modeled, 
which is expected from linear response theory. It is also 
interesting to note that in the overdamped limit for the 
monomer motion, and the continuum limit for the fluid, 
our apmpach is identical to the Oono-Freed equations of 
motionEj, which are commonly used in polymer solution 
theory. 

The main justification of our approach relies on the 
fact that a hydrodynamic (Navier-Stokes) description of 
the fluid works down to very short (actually, surprisingly 
short) length and time scales. Therefore, one should ex- 
pect that the flow around a monomer should be describ- 
able by the solution of the Navier-Stokes equation as 
soon as the distance is larger than a few lattice spacings. 
The same argument holds for the analogous MD system, 
where one expects Navier-Stokes behavior beyond a few 
particle diameters. Therefore, we may say that any two 
local couplings (for example, our LBM friction ansatz vs. 
MD) are equivalent as soon as they produce the same 
long-range flow fleld. If this is the case, then the hydro- 
dynamic interaction between two monomers (as long as 
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they are not too close) will be identical, and the single- 
monomer mobilities will also match (note that for a par- 
ticle which is pulled through the fluid at constant velocity 
by a constant force, the friction coefficient is determined 
by the energy dissipated in the surrounding flow field). 

This latter property actually allows for an easy deter- 
mination of the simulation parameter which we will 
now, for the sake of clarity, denote by the symbol Cbarc- 
A heuristic procedure, which was followed in Ref. is 
to vary this parameter in a set of simulations of a sin- 
gle monomer in solvent (which can be done very easily), 
and to measure the momomer diffusion coefficient Dq, 
until the latter has the desired value. If viscosity and 
fluid density match as well, then the long-range parts of 
the flow fields (beyond a few lattice spacings) must look 
the same. It should be noted that the Einstein relation 
Dq — ksT / C,cii thus defines an effective or renormalized 
friction coefhcient, which differs from the original bare 
one, as it contains all the backflow effects. Since these 
tend to increase the mobility, one has Ccff < Cbarc- More 
quantitatively, one can argue as follows: Let us consider 
a particle which is pulled through the solvent at constant 
velocity V by an external force F. Then, rewriting Eq. 
^, we find 



V = 



1 



Cba 



(17) 



where Uav is the flow velocity averaged over the nearest 
lattice sites of the particle, as implemented by our inter- 
polation procedure. However, to a good approximation, 
the flow field should be given by the Oseen tensor: 



(1 + ri 



(18) 



where r is the distance from the particle. Hence the aver- 
aged flow field should — in our case of averaging roughly 
at a distance a from the particle — have the form 



g-qa 



(19) 



where g is an unknown numeric constant describing the 
details of the lattice geometry and of the averaging pro- 
cedure. For example, doing the average over a sphere 
of radius d, one would directly obtain Uav — F/(67r77d), 
from which one easily derives Stokes' law. Combining 
these results and using CoffV = F one obtains 



1 

Ccff 



1 



Cba 



1 



(20) 



i. e. the overall mobility is simply the sum of the bare 
mobility and a hydrodynamic, Stokes-type contribution, 
where the lattice discretization serves to provide an ef- 
fective Stokes radius of the monomers. This relation has 
been tested by running several simulations at different 
bare couplings and different lattice constants; the agree- 
ment is remarkable, as seen from Fig. ^, where we plot 
Tya/Ccff as a function of rya/Cbarc- The parameter g is thus 
found to have the value g « 25 for our method. 




FIG. 2. Test of the predicted relation between bare and ef- 
fective friction coefficient, Eq. |2^. Grids of different lattice 
spacings were used as indicated in the figure. 

The lattice constant a hence appears not only as a pa- 
rameter which controls how accurately the Navier-Stokes 
equation is solved (this is the usual case for Navier-Stokes 
equation solvers), but it is being assigned an additional 
meaning as an effective Stokes radius. For that reason, 
it cannot be varied arbitrarily, but only within limits: A 
too small lattice constant would result in an unphysically 
large particle mobility, even if Cbare is very large. This is 
quite different from conventional Navier-Stokes equation 
solving, where one obtains systematically better results 
when a is decreased, and can be viewed as the price which 
has to be paid for introducing the simple and computa- 
tionally fast concept of a point particle, which is however, 
strictly spoken, unphysical. It should be noted that Cbarc 
controls the degree of coupling to the flow fleld: For small 
Cbarc, one has Ccff ~ Cbarc, whilc for large Cbarc the Stokes 
contribution prevails, Ccff ~ 5"7a- It should have become 
clear that hence Cbarc has no real physical meaning what- 
soever; it is really the effective friction which matters for 
the coupling. 



III. SINGLE CHAIN SIMULATION 



A. Input Parameters 



The present model is intended to represent the same. 
physical situation as an existing pure MD simulationQ. 
We therefore choose the physical input values for the new 
method as obtained by the former (all values are given in 

IIB| ). The fluid is char- 



the unit system specified in Sec. 

acterized by the the temperature ksT = 1.2, the density 
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p = 0.864, and the kinematic viscosity v = 2.8. The 
parameter /i (the fluid particle mass) is unimportant; its 
value can be absorbed in a re-definition of the n^. The 
lattice constant a of the grid is set to unity; this is roughly 
the same as the bond length of the polymer chain, and 
the interparticlc distance of the MD fluid. As in the pure 
MD simulation, we study chains of length A^ch — 30, 40 
and 60. The corresponding grid sizes (which are impor- 
tant parameters, since they determine the hydrodynamic 
interaction of the chain with its periodic images, see Ref. 
^) are L = 18, 18, and 22, respectively, which is roughly 
identical to the corresponding MD box sizes. 

The parameters for the FENE potential are taken from 
the MD simulati on as i?o = 2.0 and k = 7.0. As already 
discussed in Sec. 11 C , this does however not assure that 
the static conformations are identical: In the MD case, 
there is also the influence of the solvent, which is absent 
in the present method. Actually, the data show a sys- 
tematic deviation, which is however not very large (see 
Sec. pIBl). 



The mass of the monomers was set to unity. This ac- 
tually differs from the MD case where the monomer mass 
had been set to two. However, we also used a monomer 
of mass one in order to determine the "bare" friction co- 
efficient Cbarc , using the procedure outlined at the end 
of Sec. |ll C| , such that we found Cbarc = 20.8 from the 
requirement that the monomer diffusion coefficient has 
the value known from MD, Dq = 0.076. Had we used a 
monomer of different mass, we would also have obtained 
a slightly different value for Cbaro (these are very small 
effects, beyond what the simple picture which underlies 
Eq. 20 can capture). Since however on the time scale of 
Brownian motion it is only the parameter Coff = / Dq 
which matters, we expect an influence of the mass param- 
eter only for short times, where the dynamics differs from 
MD behavior anyways. 



It remains to specify the time steps At aad r. A choice 
of At = 0.01 is optimal for the MD partfl. Concerning 
the LBM time step t it is desirable to make it as large 
as possible because the fluid calculation is the CPU in- 
tensive part of the method. Test simulations showed the 
limiting factor to be that rii is getting negative for too 
large time steps due to increasing fluctuations, in par- 
ticular near the monomers. This situation, however, can 
always happen, although with decreasing probability for 
smaller time steps. We found that using a time step of 
T = 0.05 only approximately each lO'^th random num- 
ber one rii became negative, while for t — 0.01 such 
a case never occurred during the observation time. We 
decided to generate new random numbers in such rare 
cases, which of course slightly changes the distribution 
of the simulated noise, but is justified if the probability 
for negative n^'s is low enough. We ran the simulations 
at T = 0.05 and also did a simulation for the smallest 
system (TVch = 30, L = 18) using r — 0.01 in order to 
check the results. 

Furthermore, we should comment in some more detail 
on the lattice constant a. The choice a = 1 seems intu- 
itively reasonable, since this matches the bond length and 
the interparticlc distance in the MD system. However, 
one would in principle like to make the lattice spacing as 
large as possible, since, for constant overall volume, the 
computational effort scales as . For this reason, we 
also did a test run with a — 2 for the A'ch = 30 system, 
where we of course had re-adjusted the bare friction, see 
end of Sec. 11 C . It turned out that the decay of the 
dynamic structure factor looks quite similar. However, 



there are systematic discrepancies (see Sec. |VB|), such 
that the gain in speed is paid for by a certain loss in ac- 
curacy. In what follows we will always refer to the case 
a = 1, unless explicitly stated otherwise. 



Chain length 


30 


30 


40 


60 


LB time step r 


0.05 


0.01 


0.05 


0.05 


exponent v 


0.621 ±0.004 


0.620 ± 0.002 


0.637 ± 0.002 


0.637 ± 0.002 




94 ±5 


90 ±4 


134 ±4 


217 ± 10 


,i 


14.3 ±0.5 


13.9 ±0.4 


20.6 ±0.3 


33.5 ±0.9 




0.299 ±0.005 


0.300 ± 0.005 


0.261 ± 0.005 


0.215 ±0.004 




0.1512 


0.1525 


0.1179 


0.0986 


ksT 


1.139 ±0.003 


1.2056 ±0.003 


1.139 ±0.003 


1.139 ±0.003 




0.9951 ± 0.0004 


1.009 ± 0.0002 


1.0001 ± 0.0001 


1.006 ± 0.003 


51-exp.'' 


0.6415 ± 0.001 


0.6747 ± 0.001 


0.6630 ± 0.0006 


0.6704 ± 0.002 


DcM 


6.533 X 10"^ ± 1 X lO"'^ 


6.102 X 10"^ ± 1 X lO"'^ 


4.860 X 10"^ ± 2 X 10"^ 


3.387 X 10"^ ± 1 X 10"^ 


Do' 


0.081 


0.062 


0.076 


0.054 


Tz (estimate) 


365 


380 


705 


1650 



TABLE I. Single chain properties 



^no error due to complicated calculation 

''exponent obtained by fitting a power law in the sub-diffusive scaling regime t £ [20 ; 80] 
'^calculated using Eq. ^ 
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An important point concerning the comparison with 
analytical theory should be mentioned here. It is usually 
assumed in these theories that the time scale for the evo- 
lution of the hydrodynamic interaction is much smaller 
than the diffusion time scale of a monomer, i. e. the 
Schmidt number Sc — ^ ^ 1. This parameter can be 
set arbitrarily in our method: v is an input parameter 
and Do can be tuned by choosing Cbare- In our case, we 
have 5*0 w 32. 



B. Chain Statics 

The results for the chain lengths of A'ch — 30, 40 and 
60 are listed in Table |. The measurement of the chain's 
temperature provides a first consistency check of the al- 
gorithm. The values for fesTmoasurcd = all;;;- -^^kin show a 
discretization error of 5% for the large time step r = 0.05. 
For the small time step r = 0.01 the error decreases sig- 
nificantly. 

The radius of gyration 



1 



(21) 



with r,- 



r, |, and the end-to-end distance 



without surrounding LBM fluid. The conformations must 
be the same, i. e. the structure factors must coincide (up 
to discretization errors, which may look somewhat differ- 
ent for the chain coupled to the LBM fluid). As is seen 
from the figure, the agreement is very good, i. e. the 
method is validated to produce correct static conforma- 
tions. 



10' r ' ' ' ' ' ' " 

" N,,=30, x=0.05 



10' 



10 



10"' 



i N„,=40, x=0.05 
» N„„=60, x=0.05 
N h=30, no solvent 



10" 



10" 
k 



10' 



(22) 



are related to the number of monomers by the static ex- 
ponent ly, 



(23) 



for a self-avoiding walk ly « 0.588 from renormalization 
group theory methods and Monte Carlo simulatiorH. In 
principle, v can be obtained from the scaling law (23); 
however, this would require simulations covering a wide 
range of iVch. Hence, it is advantageous to use the static 
structure factor 



(24) 



sm{krij) 



kr- 



which probes different length scales even for a single poly- 
mer. In the scaling regime Rj^ <^ k a^^ (oq being a 
microscopic length of the order of the bond length) the 
relation 



S{k) oc fc^i/'' 



(25) 



holdsi. By fitting a power law to our data (see Fig. H) we 
get the values for v of Table | which are about 6% higher 
than the asymptotically correct value, resulting from the 
finite chain length. In Fig. ^ we also include data which 
have been generated from a simulation of a single chain 



FIG. 3. The static structure factor of the chains. 



The hydrodynamic radius 



1 



N. 



ch 



E 



(26) 



is an interesting quantity because the Kirkwood pred jc- 
tion for the diffusion of the chain's center of massESIlil 



D, 



CM 



kBT_ 

iVeh 67r77 \Rh/ , 



(27) 



depends on it. This formula, however, is only correct for 
a single chain in an infinite medium. In a finite box one 
has to take into account the hydrodynamic interaction 
with the periodic images. This will eventually lead to a 
finite-size corrected hydrodynamic radius. Quite gener- 
ally, one must expect a finite size effect of order for 
every dynamic quantity, corresponding to the slow 
decay of hydrodynamic interactions. A detailed descrip- 
tion can be found in Refs. 0,^, so that we can restrict 
ourselves to the essential points. Within the Oseen ap- 
proximation, the diffusion tensor is given by 



D, 



D(I■^,) 



k ® k 



k#0 



fc2 



exp(zk-ry) (28) 



for i j, where k = 27rn/L (n being a vector of inte- 
gers) runs over the reciprocal lattice vectors and k is a 
unit vector in the direction of k. For i — j, one has the 
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monomeric diffusion coefficient Dq, plus the contribution 
due to the hydrodynamic interaction of that bead with 
its own periodic images, 

D„; = Dol + hm ( D(r) - + f ® f) ) . (29) 

The last two expressions can be calculated efficiently us- 
ing the Ewald summation technique. The center of mass 
diffusion constant is given by 



Dc 



1 



1. 



M 



Inserting Eqs. E 

n - ^° 

JVch 



^ V -Tr(Di,, 

I and ^ one obtainsB 
2.837 keT 1 



ch2 



E 



(30) 



Tr(D,,), (31) 



which defines, by comparison with the Kirkwood formula 
(|^), a finite size corrected hydrodynamic radius: 



CM,L 



Do 



ksT / I \ 
G-Kri \Rh / 



(32) 



Rh is thus effectively increased by the periodic images. 
For our box sizes, the discrepancy between {RJj'^ ^ and 
{R^) amounts to approximately a factor of two (cf. 
Table This is in agreement with the corrections found 
in Ref. ^ 



100 



75 



^„ 50 



N^,=30, T=0.05 
N„,=40, T=0.05 
N„,=60, T=0.05 
N^,=30, x=0.01 




1000 
t (LJ units) 



2000 



The Zimm time tz, i- e. the longest relaxation time of 
the chain, is given by the condition that the chain has 
moved its own size during tz, or Dqmtz oc implying 
Tz oc which defines the dynamic exponent 2 = 3. 
This exponent then quite generally relates times to cor- 
responding lengths, such that, for example, the mean 
square displacement of a monomer on time scales below 
Tz, but above the microscopic time scales tq, should be 
proportional to t^/^ = t^/^. For a chain without hydro- 
dynamic interaction (Rouse model), where Dcm k 
one finds z = 2 + \/ v from analogous considerations. 

Figure ^ shows the mean square displacement of the 
chain's center of mass 



ch ' 



53 (t) = ( (RcM (io + t) - RcM {to)y 



(33) 



By fitting a power law we obtain the exponents and the 
diffusion constants shown in Table |. Obviously, the ex- 
ponents support the prediction of simple diffusive behav- 
ior [t^). One would expect theoretically that two diffu- 
sive regimes exist, both exhibiting behavior but dif- 
ferent prefactors, with a smooth crossover around the 
Zimm time. The accuracy of the data does not allow 
to support this crossover, which is not surprising as the 
short-time and long-time diffusion-eoastant are expected 
to be rather close to each otherll3llll'c3. In principle, the 
scaling behavior of -Dcm provides a test of the Zimm pre- 
diction DcM CK ^cJi'^- But there are large corrections-to 
scaling due to finite chain length and bead size effectsQcj. 
Therefore it is more useful to analyze the non-asymptotic 
relation (|3l| ) by comparing the values for Dq that can be 
obtained from Eq. where finite chain length and finite 
box size are taken into account, with the input value of 
= 0.076. The values are also listed in Table | show- 
ing quite reasonable agreement. Without the finite size 
corrections, the agreement is unacceptable, such that a 
negative value for Dq would be obtained. 



10" r 

. Nch=30. 1=0.05 
^ N^,=40, x=0.05 
N„,=60, t=0.05 
10 - N,=30,T=0.01 



10 



10" 



FIG. 4. The mean square displacement of the chain's cen- 
ter of mass. 



10"' 



C. Chain Dynamics 



^ -1 1 2 3 4 

10 10 10 10 10 10 

Time t (LJ units) 



The dynamic scaling picture for Zimm dynamics! FIG. 5. The mean square displacement of the central 
starts from the prediction Dcm oc R~^ (cf. Eq. |27|). monomer. 
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The mean square displacement of a single monomer i 
(which should only be evaluated for monomers near the 
center of the chain to eliminate end effects) 

5i(t) = ((r,(t + io)-r,(io))') (34) 

is plotted in Fig. ^ In the time regime below the Zimm 
time and above the ballistic regime, the scaling behav- 
ior gi{t) oc t^/^ is predicted. The corresponding fit to 
our data yields the exponents of Table |. The values ob- 
viously favor the Zimm model compared to the Rouse 
model, which predicts gi{t) oc t^/^ = t^-^^. 

The Zimm time can be estimated from the mean square 
displacement of a monomer in the center of mass system, 



52 (t) = {{[r^{t + to) - RcM(i + to)] 



[n{to)~'RcM{to)]y 



(35) 



which is depicted in Fig. y. Theoretically, a crossover 
to a plateau should evolve at the Zimm time. However, 
the crossover is quite extended in our simulation, making 
it difficult to extract a specific time for it. We therefore 
estimate the Zimm time from 



TZ 



which yields the values shown in Table 



(36) 
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FIG. 6. The mean square displacement of the central 
monomer in the chain's center of mass system. 
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FIG. 7. Normalized autocorrelation function of the Rouse 
mode Xp for different p, for the longest simulated chain 
A^ch = 60. The upper part of the figure uses Vpt as scaling 
argument, where Fp was calculated directly from the chain 
conformations. The middle part uses p^'^t, where naive dy- 
namic scaling has been applied, while the lower part also takes 
the correction factor r(p) (see Appendix^ into account. The 
meaning of symbols is the same for all three parts (see middle 
part). 



It is interesting to perform a Rouse mode ajSalysis. For 
this purpose one defines the Rouse modes aala 



JVch 

= ^ch^ cos 
n=l 



(37) 



It is well known that these modes are the (independent) 
eigenmodes of the random walk Rouse modelB. 



However, for reasons of translational symmetry along 
the chain, one must expect that the cross correlation 
(Xp(i -I- to)Xg(to)} {p 7^ l) is at any rate quite weak, 
regardless of chain statistics and dynamics, such that 
the modes can be viewed as independent modes even be- 
yond the random walk Rouse case. For a ring polymer, 
this can be shown rigorously, since in this case there is 
strict invariance under the transformation n — > n + 1, 
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such that the Rouse modes (which are then defined with 
an exp (ipTrn/Nch) factor) are eigenfunctions under this 
transformation. Hence, if end effects are not too strong, 
one should also expect for our case an independence of 
the Rouse modes. Indeed, within the accuracy of our 
data, the cross correlation terms are zero. 

Furthermore, within the approximations of the Zimm 
model, the autocorrjelation function of the modes should 
decay exponentiallyEl, 



(Xp(t + to)Xpfa)) 



exp(-t/Tp). 



(38) 



In Fig. 0, we therefore plot, for p > 1, the normalized 
autocorrelation function semi-logarithmically as a func- 
tion of properly scaled time. Firstly, we estimate Tp via 
the initial decay rate 



" " dt 



(Xp(t)X^(O)) 



(39) 



which can, within the framework of Kirkwood-Zimm the- 
ory, be calculated in terms of purely static averages, 
i. e. from the chain conformations in combination with 
a model diffusion tensor, for which we use Eqs. ^ and 
p9| . The details of this approach are described in Ap- 
pendix ^ Interestingly, it turns out that this quantity 
is only subject to an L^'^ finite size effect (which we 
neglect), in contrast to the usual L^^ behavior. This 
result holds beyond the various approximations of Ap- 
pendix our interpretation is that any contribution of 
global center-of-mass motion of the chain is being sub- 
tracted, such that the leading-order hydrodynamic in- 
teraction with the periodic images cancels out, and only 
a dipole-type interaction remains. In the upper part of 
Fig. ^ we thus plot the autocorrelation as a function of 
Tpt, where Tp was calculated directly from the simulated 
chain conformations, in combination with the Oseen ten- 
sor. It is seen that the Oseen formula describes the decay 
quite well; however, the data collapse is not particularly 
good. There is also some curvature, indicating a non- 
exponential decay. The middle part of the figure then 
uses the scaling argument p^'^t. This p-dependence re- 
sults from the calculation of Tp, where instead of the ac- 
tual chain conformations asymptotic self-avoiding walk 
statistics is employed (see Appendix ^), as the leading 
power law. This corresponds to simple dynamic scaling, 
which views the pth mode as equivalent to a chain of 
length Nch/p, such that Tp oc [Nch/vY^ ■ However, the 
more detailed calculation of Appendix ^ yields an addi- 
tional weak p-dependence, i. e. a correction factor r(p), 
whose presence indicates, in our opinion, that the simple 
picture of subchains of length N^h/p is not fully justified. 
Taking this correction into account, we obtain a very 
nice data collapse (see lower part of Fig. This is quite 
remarkable; one would of course expect the best data col- 
lapse for the uppermost part which involves the smallest 
number of approximations. It seems that there are vari- 
ous errors involved which somehow happen to cancel out. 
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FIG. 8. Scaling plot of the dynamic structure factor for 
A^ch = 60, for Rouse scaling (z = 3.7), asymptotic Zimm scal- 
ing (z = 3), and z = 2.8, which produces the best collapse. 

As far as the absolute value of the decay rate is 
concerned we find reasonable agreement: While the 
lower part of Fig. ^ shows a decay rate of roughly 



3 X 10 ^p'^''r(p), Eq. |A21 predicts a decay rate of or- 



der 5.4 X 10~ p r(p) where we have for simplicity used 
the random walk value for the constant A, and = 2.0 
(extracted from the results for via — b'^N^'^). 
The dynamic structure factor 

= 1^ E (e^P ■ [^^(*) - (0)])) (40) 



10 



is predictedi to exhibit the scahng behavior 
S(k,t) = S{k,0)f{k't) 



(41) 



if both wavenumber and time are in the scahng regime, 
i. e. ^ fc ^ flg ^ and tq <^ t <^ tz- It is even 

possible to calculate explicit formulas (rigorously for the 
random walk Rouse model and usiiici-4he linearization 
approximation in the Zimm case)l3ca'E3, which suggest 
that there is an exponential dependency on (fc^i)^/^ for 
Tkt ^ 1, where Tk is the (fc-dependent) decay rate. 
Hence a plot of S{k,t)k^/'' against (k^t)"^/^ should — 
for the correct model — collapse to a straight line in a 
log-linear representation. For iVch = 60, the results are 
shown in Fig. || (the plots for the other chain lengths look 
quite similar). The data were restricted to the scaling 
regime 20 < t < 80 and 0.7 < fc < 2. These ranges were 
obtained from the single-monomer mean square displace- 
ment. Fig. |[ and from the static structure factor. Fig. ||, 
respectively. Values of S{k,t) below 0.01 were discarded, 
for reasons of statistical accuracy. It is clearly visible 
that the simulation shows Zimm rather than Rouse be- 
havior. A dynamic exponent of z = 2.8 yields the best 
data collapse. Such an effective value, which is, due to 
corrections to scaling, somewhat smaller than the cor- 
rect asymptotic one, is quite usually observed, not only 
in simulationsO, but also in experimentsa. 

Concerning finite size effects, one has for a finite box 
size S = S{k, t, L), and scahng is corrupted by the second 
length L in the problem. The influence can be estimated, 
in close analogy to the procedure presented in Appendix 
by studying thep4kGasu formula for the fc-dependent 
diffusion coefRciente3'cZl, 



D{k,L) 



Ey (k-Dij •kexp(ikry ) 

(exp(ikrij)) 



(42) 



which is L-dependent because of the finite size form ( p^ ) 
of Dij. D{k,L) is related to the initial slope of the dy- 
namic structure factor via 



D{k,L) 



1 



f S{k,t,L) 



J™ kH^''' yS{k,0,L) 



(43) 



We do not present the details of our semi-quantitative 
analysis here since they have been outlined in Ref. |7| al- 
ready. The result is a fc-independent correction term 
of order (note that S{k,t) does contain the overall 
chain motion, for every wavenumber). As the leading- 
order (L = oo) term is proportional to k in the scaling 
regime, the conclusion is that scaling is corrupted, but 
the relative contribution of the finite size correction gets 
weaker with increasing k. For the k —>■ limit finite size 
corrections amount to roughly 100%, as has been shown 
by the calculations for the hydrodynamic radius in Sec. 
[II B. In the scaling regime the corrections are much 
smaller, because it is closer to the kL — > oo limit. This 
is, in our opinion, the main reason why the data collapse 
works so nicely. 



IV. COMPARISON TO THE CORRESPONDING 
MD SYSTEM 



A. Efficiency 



Since the system is highly dilute, the CPU cost for the 
MD part for the polymer chain is negligible, and the lat- 
tice Boltzmann part uses up practically all computational 
resources. It should be noted that this part can be op- 
timized by choosing appropriate simulation parameters; 
our choice {a = 1, t = 0.05) is probably not the most 
efficient one. Firstly, it is possible to increase the lattice 
spacing somewhat, without substantial loss in accuracy. 
For example, going from a = 1 to a = 2 reduces the com- 
putational effort by a factor of eight. This increase seems 
howe ver t o be slightly too large already; as outlined in 
Sec. IVB, a = 2 produces less accurate data. Secondly, 
one can try to exploit Eq. ^ by varying r or a while keep- 
ing 1/ — 2.8, such that the simulation runs at A = —1, 



for which the LBM algorithm takes a particularly simple 
form in. which a substantial number of operations can be 
savedE3. Further speedup can be expected if the require- 
ment 1/ = 2.8 and Dq = 0.076 (for mapping to MD) were 
released. However, we have not checked these questions 
in a systematic fashion; in particular, our discussion has 
not taken into account that the limit of stable time steps 
T depends on both a and A in a non-trivial way. We 
hence want to simply state that our present choice of pa- 
rameters is not yet a fully optimized one; therefore the 
numbers given below (for a = 1 and t = 0.05) should 
be viewed as a lower bound of the efficiency which the 
method can attain. 

On one EV5.6 processor of a 433 MHz DEC Alpha 
server 8400 (for a typical box size of L = 40) our code 
obtains 3.1 x 10^ grid point updates per second. In 
order to compare this number with the molecular dy- 
namics system, we note that one grid point corresponds 
to 0.86 solvent particles for p = 0.86 and a — 1.0. 
Therefore, the efficiency of the code in MD units is 
3.1 X 10^ X 0.86 « 2.7 x 10^ particle updates per second. 
This number should be contrasted with the efficiency pf 
optimized MD codes for short-range LJ fluids, which isc£l 
(on the same machine) 2.1 x 10^ particle updates per sec- 
ond, using the code described in Ref. 49, Thus, the LBM 



would run by a factor of 1.3 faster than MD if the same 
time step were used. However, the lattice Boltzmann 
time step r — 0.05 is more than an order of magnitude 
larger than for the pure MD system: The latter must be 
run without friction and noise, i. e. in the microcanonical 
ensemble, in order to strictly conserve momentum (other- 
wise the hydrodynamic interaction would be screenedEJ) . 
Such a simulation can only be stable on long time scales 
if the timej-step is sufficiently small; according to our 
experienceaU, one needs At = 0.003. Taking these fac- 
tors into account, we obtain a net speedup of a factor of 
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22, which, as outlined above, can be increased further by 
choosing a coarser lattice, i. e. by trading in accuracy 
for speed- A detailed comparison with the "competitor" 
DPDEZra is highly desirable, but not done here, last not 
least becausc-the match of the viscosity is much less triv- 
ial in DPDcJO. From what we know from the literature, 
we expect that the two methods would be roughly com- 
parable in speed, at least by order of magnitude. 
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FIG. 9. The dynamic structure factor S{k, t) for the new 
method with r = 0.05 (circles) and r = 0.01 (line) compared 
to pure MD simulation (crosses) for three different k values 
(iVch = 30). 



B. Static and Dynamic Behavior 

In order to check how well the new method produces 
the same physics as the original MD modefl, from which 
all simulation parameters were derived, we focus on the 



comparison of the structure factor S{k,t) for both meth- 
ods, as shown in Fig. jo] (time dependence at constant kY 
Fig. |l^ (fc dependence at constant time), and Fig. |ll| 
(time dependence for the normalized structure factor). 

Let us first consider the static case t — 0. The corre- 
sponding plot (Fig. ^ uppermost part) for A^ch = 30 
shows systematic deviations. These manifest for exam- 
ple in the discrepancies of the static scaling exponent 
{i^ — 0.59 for the pure MD simulation, v = 0.62 . . . 0.64 
for the new method); the chain is more stretched using 
the new method. The absolute values for the static struc- 
ture factor differ up to about 25%. Similar results hold 
if one compares other static quantities like the radius of 
gyration or the end-to-end distance. It can be verified 
that the discrepancies show no significant dependency 
on the chain length for the range investigated here (30- 
60). Moreover, they are not due to a discretization error 
in time, as the plots for t = 0.05 and r = 0.01 show. 
The reason is rather simply the fact that the MD chain 
is subject to a different potential (intra-chain plus sol- 
vent) than the LBM chain (intra-chain only). For that 
reason, there is a systematic difference in the static con- 
formations, which then, in turn, will also affect the dy- 
namic properties somewhat. For example, the Nch = 60 

1/2 

chain has a gyration radius — 5.79, while the 

corresponding MD chainQ has a gyration radius of only 
4.78. It is hence not surprising that the larger chain is 
also somewhat slower, as the comparison of the diffusion 
constants confirms {Dcm — 3.39 x lO^'^ for the larger 
LBM chain, and Dcm = 4.25 x 10"^ for the smaller MD 
chain). Therefore, in order to achieve a better match of 
static and dynamic properties, it would be necessary to 
re-adjust the potential for the LBM chain such that the 
conformations are more similar. This is possible, but not 
completely trivial, and has not been attempted in this 
work. On the other hand, for the dynamics parameters 
(i. e. the viscosity and the friction coefficient), it is quite 
easy to achieve matching, as has been described in Sec. 
[II A. 



Turning to the decay of S{k,t), we first note that the 
direct comparison of the data (Fig. || and |l^) yields sim- 
ilar discrepancies of up to 25% as for the static case. The 
overall agreement is however quite reasonable. In order 
to divide out the trivial amplitude effect, we also plot 
S{k, t) /S{k, 0) for three different k values in Fig. ^ For 
k in the scaling regime, the agreement is much better, 
with differences of a few percent only. This is not too 
surprising, since in this regime the decay rate should in 
essence be given by k^ksT/r] times a numerical prefactor 
which depends only weakly on the details of the chain 
statisticafj'cZl. In the long-wavelength regime (inset of 
Fig. |ll|) the decay is given by exp {—Dcuk^t), which is 
nicely confirmed by the data, and thus the ratio of the 
decay rates is just the ratio of the diffusion constants, 
i. e. there is again a discrepancy of roughly 20 % (this is 
hardly visible in Fig. |l^, due to noise in the MD data). 
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FIG. 10. The dynamic structure factor S{k, t) for the new 
method with r = 0.05 (circles) and r = 0.01 (hne) com- 
pared to pure MD simulation (crosses) for three different 
times (iVch = 30). 

To summarize, we find that both methods are well- 
suited for quantitatively reproducing the dynamics of 
polymer chains in solvent, and both reveal Zimm behav- 
ior very nicely. The discrepancies which we find in the 
dynamic properties can be directly traced back to the 
non-perfect match of the static conformations. If those 
had been matched by an adjustment of the potential, 
then the agreement would probably be close to perfect. 

Finally, let us discuss in more quantitative terms the 
influence of the lattice spacing. To this end. Fig. |l^ 
compares the decay of the normalized dynamic structure 
factor of an A'ch = 30 chain for three k values, obtained 
by running the same system with two different lattice 



spacings a = 1 (as discussed previously) and a = 2. All 
other simulation parameters (in particular the box vol- 
ume, and the monomeric diffusion coefficient Dq — not 
the bare coupling Cbarc) were left identical. As one sees 
from the figure, the larger lattice spacing induces decays 
which are systematically slower, by roughly 20 % to 25%. 
It is thus a question of desired accuracy if one wants to 
consider these results as still acceptable or not. The ob- 
served effect goes in the direction which one expects, for 
the following reasons: As soon as the lattice spacing ex- 
ceeds the size of the chain, there will be no hydrodynam- 
ics left and one will observe pure Rouse dynamics, which 
is slower. Of course, this must be a systematic crossover 
as a function of lattice spacing. Thus one expects a de- 
crease of the hydrodynamic correlations with increasing 
a (also consistent with the reasoning at the end of Sec. 



II C ), and hence a systematic slowdown of the dynamics. 
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FIG. 11. S{k, i) /S{k, 0) for A^ch = 60 using the new method 
(solid lines) with t = 0.05 compared to pure MD simulation 
(dashed lines) for three different k values. 



V. CONCLUSION AND OUTLOOK 

With this paper we have established a new method to 
simulate polymer-solvent systems. The solvent is mod- 
eled by the lattice Boltzmann method and the polymer 
by a continuum bead-spring model. The two parts are 
coupled using a simple dissipative friction ansatz which 
locally conserves mass and momentum. The driving force 
of the system are thermal fluctuations which are added 
to both the fluid and the polymer. The main advantage 
of the new method compared to MD is its computational 
efficiency, which amounts to a factor of 20, or even more, 
if one is willing to be s atisfie d with less accurate results. 

As described in Sec. Ill A, it is possible to obtain the 
physical input parameters for the new method from re- 
sults of existing MD simulations. Therefore, one can view 
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the present method as a coarse-graining procedure where 
one goes in a well-controlled way from small length and 
time scales to larger ones. As the results show, this is 
possible without substantial loss of information about the 
statics and dynamics on the mesoscopic scale. 

The input which is needed from a more microscopic ap- 
proach consists of: (i) Effective potentials for the coarse- 
grained monomers such that the static chain conforma- 
tions are roughly reproduced (this was the part to which 
we did not pay much attention, with the result that this 
is the largest source for the observed deviations); (ii) the 
solvent temperature, density and viscosity, and (iii) the 
monomeric diffusion coefficient, from which one adjusts 
the coupling. 

It seems that a lattice spacing which roughly matches 
the chain's bond length and the interparticle distance of 
the solvent is optimal. A lattice constant which is chosen 
too large will result in underestimated hydrodynamic in- 
teractions, as seen from the data with a = 2, while a too 
small lattice spacing will result in a large computational 
effort, plus (if it becomes very small) a monomeric diffu- 
sion coefficient which will exceed any realistic value, due 
to an effective Stokes radius which is too small. 

We have chosen the parameters of Ref. 0for our simula- 
tion and performed a detailed quantitative comparison of 
the results. The main deviations result from insufficient 
match of the static conformations. The current model 
is therefore as appropriate as the original MD model for 
verification of Zimm dynamics in dilute polymer solu- 
tions. The dynamic scaling laws (in particular the k^t 
decay of the dynamic structure factor) could be observed, 
and there is good agreement with the decay rates pre- 
dicted by the Zimm model, if the finite box size effects 
are taken into account. Interestingly enough, the decay 
of the Rouse modes is only subject to an finite size 
effect, while most other decay rates have a large L"-'^ fi- 
nite size correction, due to the behavior of the Oseen 
tensor. 

After having tested the method successfully future 
work can now deal with more controversial problems, 
like the influence of hydrodynamics on the motion of 
a semi-flexible chain or the hydrodynamic screening in 
semidliute solutions. It should however be kept in mind 
that the algorithm in its current version is only suit- 
able for problems where the polymer concentration is 
low. The coupling only takes into account the momen- 
tum transfer between monomers and solvent. Excluded- 
volume effects between solvent particles and monomers, 
which are very important for processes like, e. g., the pen- 
etration of solvent into a dense polymer matrix, are not 
properly modeled. A study of such topics would require 
a generalization of the algorithm which would assign a 
flnite volume to the monomers. 

It is a pleasure to thank Ralf Everaers and Alexander 
Kolb for helpful discussions, and the latter for a critical 
reading of the manuscript. 
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FIG. 12. Comparison of S{k,t)/ S{k,0) for N^i, = 30, two 
different lattice spacings, and three k values as indicated. 



APPENDIX A: INITIAL DECAY RATE OF 
ROUSE MODES 

In this appendix we outline the details of the calcu- 
lation of Fp, i. e. the initial decay rate of the autocor- 
relation function of the Rouse modes for p > 1, where 
we treat the general case of a chain whose statistics is 
described by an exponent v (i. e. v = 0.5 for a random 
walk (RW), and = 0.6 for a self-avoiding walk (SAW)|)|. 
We start by stating the result of linear response theoryO, 



£ / (Xp(t)Xp(O)) 
dt I (X^) 



(Al) 



1 



E 



(X2) . \ dr, 



dr 



i/3 



where Greek indices again denote Cartesian coordinates. 
Evaluating the derivatives of the Rouse modes, one ob- 
tains 



1 ycos('^(z-l/2) 



(A2) 



cos(£-(j-l/2))Tr(D,;,} 



From the deflnition of the Rouse modes, Eq. 37 
flnds 



-rjKos I ^(i - 1/2) 

ch— V^ch 



(A3) 



cos, -(,-1/2)1, 
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which is evaluated via (6 is the bond length) 



(A4) 
(A5) 
(A6) 



(note that the last relation holds only asymptotically for 
large \i — j\)- Approximation by an integral yields 



in) 



2^ch ^0 



dx / dy\x-~ y\ 

JQ 

m pns 



2iy 



(A7) 



Furthermore, we use the relation 

cos a cos /3 = 2 [*^°^ ~ /3) + cos (a + /?)] 
and transform to the variables 

u=^[x-y), v^—[x + y). 



(A8) 



(A9) 



Exploiting the symmetry of the integrand with respect 
to u, and performing the integration over v, we find 



f{p) 



(AlO) 



with 



f{p) — — / duu^'' [sinu — (pTT — u) cosu] . (All) 

PTT Jo 

For the RW case, f{p) is exactly unity, while for the SAW 
case a weak dependence on p remains; however, also in 
this case f{p) is close to one. Using the MAPLE software 
package, we have numerically evaluated this function; for 
the first 20 Rouse modes it is tabulated in Table |l[ 

The calculation of the numerator of Eq. |A2| is per- 
formed using precisely the same procedure, the only dif- 
ference being that (^{vi — fj)^^ is replaced by Tr(Dy), 
which we calculate using the finite box size form, Eq. |2^ : 

TV (Dij) = / dk (exp(ik • ry )) , (A12) 

W Jko 

where we have replaced the summation over wavenum- 
bers by an integral 



1 



k#0 



(27r)3 



(A13) 



/co = 27r/L denoting the cutoff wavenumber due to the 
finite box size. 



2i> 



The factor (exp(zk ■ Ty )) describes tl 
chain, and must, for reasons of scaling 
have the form 

(exp(ik-ry)) ^g{k%'^\i-j\ 

It should be noted that, for reasons of inflection symme- 
try, g must depend on fc^, and that g(0) = 1. We further 
introduce the constants 



structure of the 
, asymptotically 

(A14) 



A = 



B 



dwg['uP') 
dg{w'^) 



dw"^ 



(A15) 
(A16) 



u)=0 



for a 



random walk one has g = 
g(w'^) = exp(-w^/6), A — 



1. e. 



For example, 
exp(-(&V6)fc2 

y/3TT/2, B = -1/6. We now calculate Tr(Dy ) by per- 
forming a Taylor expansion with respect to /cq = 0{L^^); 
the result is 



Tr(D 



A 



b\i 



(A17) 



Interestingly, the linear term does not depend on the 
monomer indices at all. From this, we conclude that the 
linear contri butio n to the decay rate exactly van- 
ishes, due to Eq. 

" In 



A5, and that the leading order finite 
size effect is actually of order L"^, i. e. quite small. 



what follows we will therefore only concentrate on the 
leading-order term for an infinite box. Using the same 
procedure as for (X^), one finds 



p 


RW 

h{p) = r{p) 


fip) 


SAW 
h{p) 


r{p) 


1 


1.040901 


1.229939 


1.531335 


1.245049 


2 


1.155368 


1.096321 


1.671897 


1.525007 


3 


1.186325 


1.099453 


1.711021 


1.556248 


4 


1.203640 


1.075431 


1.732468 


1.610952 


5 


1.213328 


1.077224 


1.744639 


1.619569 


6 


1.220118 


1.067140 


1.753074 


1.642778 


7 


1.224789 


1.068286 


1.758929 


1.646496 


8 


1.228399 


1.062691 


1.763420 


1.659391 


9 


1.231138 


1.063494 


1.766850 


1.661363 


10 


1.233376 


1.059915 


1.769636 


1.669601 


11 


1.235174 


1.060514 


1.771886 


1.670781 


12 


1.236696 


1.058018 


1.773782 


1.676515 


13 


1.237967 


1.058484 


1.775371 


1.677278 


14 


1.239069 


1.056638 


1.776745 


1.681508 


15 


1.240014 


1.057013 


1.777926 


1.682029 


16 


1.240849 


1.055589 


1.778967 


1.685283 


17 


1.241579 


1.055899 


1.779880 


1.685654 


18 


1.242234 


1.054765 


1.780696 


1.688239 


19 


1.242815 


1.055026 


1.781422 


1.688510 


20 


1.243341 


1.054101 


1.782079 


1.690615 



TABLE IL The functions f{p), h{p), and r(p), as defined 
in the text, for both the RW and the SAW case. 
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/ dx / dy^ rjj cos —a; cos — y 



N, 



2-u 
ch 



(A18) 
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This function also exhibits a weak p-dependence, see Ta- 
ble m even for a RW. Finally, introducing 



r{p) - h{p)/f{p), 



(A20) 



also tabulated in Table ||, we can write the result for Tp 



as 



r, = A^'^(^y\{p). 



(A21) 



The leading power-law dependence on p and A'ch is ex- 
actly what one expects from dynamic scaling. The func- 
tion r(p) is a correction to scaling. As far as the numer- 
ical prefactor is concerned, we get (in the RW case) a 
relaxation which is roughly the same, as that calculated 
in the textbook by Doi and EdwardsH. 
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